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Preface 

This  study  Investigates  the  use  of  active  coju.ro 1 
for  aircraft  landing  gear  to  minimize  runway  imposed  vib¬ 
ration.  The  optimal  control  appro-  "1  }  tai?en,  wi  fh  both 
the  (1-cos  tft)  hump  and  runway  speotivl  density 
considered  for  inputs  to  a  tv:o  degree  of  freed ov  linear 
landing  gear  model.  The  conjugate  gradient  numeri  ex  i 
technique  is  used  for  problem  solution. 

I  wish  to  extend  my  slncerest  appreciation  to 
Mr.  Ronald  0.  Anderson  at..  Captain  James  Di How  of  the 
Control  Criteria  Branch,  Air  force  flight  Dynamics 
Laboratory  for  making  the  completion  of  this  study  possii 1< 
I  wish  to  thank  also  Major  John  C.  Schoco,  my  faculty 
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Finally,  I  wish  to  thank  my  wife  find  children  for 
accepting  a  part-time  husband  and  father  for  the  duration 
of  this  study. 


Ronald  A.  De  Yoe 
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At*  a  r»fr» 

- N.*X 

The  purpose  of  this  thesis  is  to  Investigate  the  use 
of  active  control  on  aircraft  lancing:  rear.  The  problem 
of  vlibratioh  isolation  Is  considered  using  the  optimal 
control  approach  to  establish  the  control  force  ned’essary 
to  minimize  runv:ay  Imposed  vibration.  A  two  degree  of 
’freedorn  linear  landing  gear  model  is  considered,  and.  both 
the  (1-cos  v:t)  bump  and  an  equivalent  deterministic  input 
derived  from  runway  spectra,  are  used  as  vertical  forcing 
function  inputs  to  the  model.  An  integral  square  cost 
function  is  formed  with  the  Integrand  consisting  of  the 
absolute  acceleration  of  the  aircraft  squared,  plus  a 

V/fil  ciltlvicr  f *5>r» +-<"»■»*  t!i«OP  4-W*o  - 1 - '  «  • 

—  _  '  W  WM—  v  WA.  VW««VA  VA  A  V/X  V/^.  UCi  1  CU  * 

The  conjugate  gradient  technique  is  used  for  numerical 
solution  of  the  problem.  > 

Results  indicate  that  considerable  improvement  in 

i  . 

vibration  isolation  could  be  gained  with  the  active 
control  landing  gear  system  over  the  present  passive 
systems  foremost  paved  run’^ay  surfaces.  It  is  recom¬ 
mended  that  further  study  be  made  with  a  nonlinear 
landing  gear  model  and  also  with  the  relative  displacement 
of  the  model  added  to  the  cost  formulation. 


cg  c/’i-:/70-6 


Active'  Control  For 
■  Aircraft  Landing' Gear 

'  I .  Introduction 

Background 

Whenever  aircraft  make  contact  with  a  runway, 

i 

whether  Vt  be  taxiing,  take-off  or  landing,  they  are 
subjected  to  many  s true ten  1  "tresses  and  strains.  These 
stresses  end  strains  on  the  air  fra-.?  are  caused  by  forces 
transmitted  throe. h  the  .shock  strut  or  landing  gear  by 
imperfections  in  the  rmv.ov;  surface.  The'  large,  heavy, 
flexible,  airplanes,  such  as  the  Cl4l,  352,  and  C-5A,  in 

particular  pre.sent  unique  problems  of  structural  fatigue 

i 

i 

because  of  their  combination  of  increased  size  and 
increased  structuul  flexibility.  When  a  352  moves  over 
a  runttny  the  wings  move  up  and  dovm  in  an  oscillatory 

motion  due  to  the  bumps  and  indentations  in  the  runway 

/'  -\  .  ’  .*  *’ 

(surface,  and/ although  less  noticeable  to  the  observer, 

/  :  • 

the,  entire  airframe  of  the  airplane  is  subjected  to  these 
same  vibrations  from  runway  roughness  or  uneveness, 

j 

•r obi cm 

kost  aircraft  landing  gear  or  shock  struts  use  what 
Is  commonly  known  as  an  "oleo"  to  Isolate  the  aircraft 
from  runway  imposed  shock  and  vibration.  The  oleo  can  be 
considered  n  passive  spring  and  dnshpot  system.  By 
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passive.  It  Is  meant  that  the  shock  strut  simply  "reacts" 

1 1 

to  anje-shock  or  vibration  Input  from  the  runway,  apd  it 

/ 

is  not  subject  to  any  outside  or  external  control/.  The 
oleo  Is  also  primarily  designed  to  absorb  landing  impact 

ene  ray  and  herein  lies  a  problem.  The  shock  from  land  in?- 

/ 

Impact  and  the  vibration  from  runway  roughness  are  two 
distinct  environments.  Landing  impact  energy  can  be  much 
greater  in  magni  tude  than  runway  brat ion  Energy ;  /thus , 
if  emphasis  in  design  of  the  oleo  is  placed  on  absorption 
of  landing  impact  energy,  the  oleo  is  then  too  hard  of  a 
suspension  system  to  efficiently  absorb  vibration  energy 
from  runway  roughness.  In  addition,  the  input  to  the 
shock  strut  from  land  In:'  impact  is  a  well  defined 

/  i 

deterministic  input,  while  the^v-hbrationai  input  from 

a  **~  -  — ' 

runway  uneveness  varies  witfidlf  ferent  Surfaces  and  is 

/  2 

/ 

random  in  nature.  An  optimum  passive  system  can  be 
designed  for  absorption  of  landing  impact  energy,  but  it 
is  more  difficult  to  design  such  a  system  to  absorb 
random  runway  vibration  energy  (Ref  3:“9). 

r 

A  number  of  studies  such  as  those  of  Ref (1-5) 
have  shown  that  considerable  improvement  in  performance 

can  be  gained  by  using  an  active  .shock  and  vibration  \ 

!  /  '  \ 

isolation  system, the  Improvement  in  performance'  being  j 

less  vibration  transmitted  to  the  isolated  body  by  the 

/ 

/ 

active  suspension  system  with  the  same  amount  of  sus¬ 
pension  clearance  or  "rattle  space"  as  that  of  the  passive 
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suspension  system.  In  terms  of  the  aircraft,  this  would 
mean  less  vibration  transmitted  to  the  airframe  with  the 
active  system  given  the  same  suspension  clearance  as 
that  of  the  passive  shock  strut  System  or  olefo. 

The  active  suspension  system  Is  a  feedback  control 


system.  With  the  active  sy stem,  sensing— device s  tiould 
be  used  to  rive  a  preview  reading  of  the  change  in 
runway  height  from  somjr  mean  value.  Other  parameters 
such  as  aircraft  sp^ed,  w«i  ,->•>,  t,  wing  lift  etc.  would 
be  then  fed  to  an  oh- board  computer.  The  computer  would 
produce  shoc’^^strut  control  signals  which  minimize  Induced 
loads  and  vibration  from  runway  uneveness. 

When  the  sensing  system  felt  or  spotted  a  depression 
m  the  runway,  the  shock  strut  would  be^extenTTecT- just  the 
right  amount  to  maintain  the^ohe^'g"  force  or  aircraft's 
wein-ht  at  all  t  ijnes^ihls  would  tend  to  hold  the  airplane 
at  its  initial  level.  Similarly  when  the  airplane  passes 
over  a  bump,  the  shock  strut  would  be  retracted  the  proper 

amount  to  again  maintain  the  one  "g"  force  and  hold  the 

/  • 

airplane  at  its  initial  level. 

Purpose  and  Scope 

The  purpose  of  this  study  is  to  investigate  the  use 
of  active  control  on  aircraft  landing  gear.  The  objective 
is  to  find  the  control  force  which  when  applied  to  the 
shock  strut  or  landing  gear,  will  minimize  runway  imposed 
vibration.  It  is  necessary  to  first  establish  the  nature 


J 
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I  s  • 

/  i 

J  of  the  control  force  in  order  to  determine  if  such  an 

i 

j  actuator  v;hether  hydraulic,  pneumatic,  or  whatever  could 
do  the  Job.  Secondly,  it  is  desired  to  investigate  the 
optimal  control  approach  to  the  problem  with  the  itegral 
square  value  of  absolute  acceleration  and  control  as  a 
cost  criterion.  The  numerical  solution  Kill  be  obtained 
using  the  conjugate  gradient  technique  developed  by 
Lasden,  Titter  and  Warren  (Ref  8).  . 


■j- 
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T  T  .  Bynnmlc  "odel 

The  dynamic  model  to  use  for  this  kind  of  study 
is  not  readily  apparent.  The  choice  of  exact  dynamics 
of  a  landing  ge^r  of  a  particular  airplane  does  not  lend 
Itself  toward  a  feasible  study  at  this  time.  This  is 
partly  due  to  the  fact  that  there  are  so,  many  undefined 
variables  such  as  what  kind  of  actuator,  how  it  is  to  be 
used  and  any  limitations  resulting  from  a  particular 
application  of  such  an  .actuator.  It  therefore  seems 
reasonable  to  keep  the  model  simple  and  general  enough 
to  be  able  to  draw  some  Meaningful  conclusions.  The 
one  reservation  that. -ust  be  kept  in  mind  is  that  if  the 
general  model  restricts  the  number  of  degrees  of  free¬ 
dom  or  makes  assumptions  about  the  environment,  then 
there  is  trie  possibility  that  the  optimum  solution  for 
the  general  model  may  be  less  than  optimum  for  a 
particular  application. 

The  active  system  chosen^to  represent  the  landing 
gear  Is  the  two  degree  of  freedom  model  shown  In  Figure  1. 
At  first,  the  study  began  with  a  model  exactly  like  that 
of  Figure  1  except  it  had  no  passive  spring  element  Ks  or 

damper  Cs.  This  meant  the  only  suspension  element  was 

! 

the  active  device  such  as  a  hydraulic  or  pneumatic 
actuator.  But  in  examining  this-klnd  of  system  more 
closely  it  can  be  seen  that  for  isolation  of  the  air¬ 
craft  from  runway  roughness,  the  actuator  would  have  to 

t 

/ 
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Dynamic  ''odel 
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supply  a  constant  "FI"  force  under  variable  displace¬ 
ment  required  by  the  runway  bumps  and  indentations. 

The  state  of  the  art  of  this  kind  of  actuator  is  not 
advanced  enough  to  do  the  job.  However,  by  adding  the 
sprint  Xs  and  damper  Cs  as  in  Figure  1,  the  problem  takes 
on  a  completely  different  form.  The  actuator  must  supply 
only  the  additional  force  required  together  t^lth  the 
spring  8nd  damper  to  maintain  the  FI  force  on  the 
airplane. 

Referring  to  Figure  1,  the  aircraft  is  the  rigid 

mass  ’1.  The  land ing  gear  suspension  is  Idealized  as  a 

massless  element  providing  forces  between  Ml  and  tihe 

unsprung  mass  V2.  The  total  suspension  force  is  then 

the  sum  of  the  forces  produced  by  the  sorinc  Ks.  the 

damper  Cs,  and  the  actuator  Fs.  I;2  is  supported  by  the 

•  • 

linear  tire  spring  "t ,  and  the  system  is  excited  by  the 
bump  F ( t ) .  The  displacement  of  Ml  measured  relative 
to  the  fixed  frame  of  reference,  is  denoted  by  Z.  The* 
displacement  of  Ml  relative  to  M2  is  denoted  by  XI,  and 
the  displacement  of  M2  relative  to  the  ground  or  fixed 
reference  is  denoted  by  X3. 

System  Equations  — 1 

The  equations  of  motion  for  the  suspension  system  are 

(1) 

Mi  I  =  Fs  -Ks  X,  -Cs  X.  -n,q 
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Ma  X3  =  -  F 5  +-  Ks  X,  +  CsX, 

(2) 

+-.K+  (  FA)  -  Xs)-M 3Cj 

with  the  Kinematic  conditions 


z 

~  X. 

Xs 

# 

z 

$ 

=  X. 

+ 

(3) 

I  » 

z 

1  4 

=  X, 

+ 

«  * 

Xa 

Assumption? 

First,  the  runway  whether  it  be  concrete  or  soil  is 
assumed  to  be  rigid,  and  the  Is  riding  gear  wheel  is 

a» 

restrained  to  follow  the  runway  profile  (  no  wheel  hop 
Is  permissible  ).  Secondly,  the  model  Is  a  rigid  mass 
model  and  will  predict  only  low  frequency  behavior,  but 
the  frequencies  of  interest  are  fairly  low  and  are  in  the 
range  of  ,5  to  32  Hz.  This  frequency  range  represents 
many  of  the  important  natural  modes  of  oscillation  of 
aircraft  such  as  the  B52,  and  C-5A.  The  third  assumption 
is  of  course  for  the  active  landing  gear  system  in  this 
study,  it  Is  assumed  that  the  characteristics  of  the  bump 
on  the  runway  are  known  before  the  landing  gear  wheel 
rides  over  it.  In  other  words,  some  kind  of  a  runway 
preview  scanning  system  is  used.  Finally,  aircraft 

landing  gear  generally  allow  for  more  than  two  degrees 

/ 

\ 

R 


of  freedom;  however  only  vertical  notion  is  considered 
since  the  model  is  Intended  primarily  for  study  of 
aircraft  response  to  vertical  runway  excitation. 

System  Inputs 

Two  types  of  input  are  used  for  F(t).  The  first 
Is  the  (1-cos  wt )  bump,  where  the  frequency  v;  is  deter¬ 
mined  by  dividing  the  velocity  of  the  aircraft  by  the 
wave  length  A  of  the  burp  and  then  multiplying  by  2 IT  . 
The  argument  of  the  cosine  function  is  less  than  or  equal 
to  2TT and  for  this  input,  the  equations  of  motion  are 
solved  for  the  ti"*~  duration  of  the  bump.  The  second 
type  of  input  Is  on  equivalent  deterministic  Input 
derived  from  rurr-'e.y  spectra.  It  is  a  special  function 
such  t-.vvai-  ty.a  ww.'^  rev.ghn?cc  ic  rsprcswitcd  =». 
random  variable  and  root  mean  square  system  output  values 
of  force,  displacement,  velocity  etc.  can  be  analyzed. 
This  input  takes  the  form  of  a  decaying  exponential  and 
is  derived  in  Appendix  B.  For  this  kind  of  input,  the 
equations  of  motion  are  solved  over  a  long  enough  time 

i  • 

period  until  no  change  in  system  response  is  detected 
by  increasing  the  time  period. 

The  decaying  exponential  derived  from  runvzay  spectra 
is  advantageous  in  that  it  does  not  require  a  range  of 
frequencies  to  be  covered  as  does  the  (1-cos  wt)  bump, 
and  it  provides  a  means  to  obtain  root  mean  square  system 
output  values.  However,  runwr-y  spectral  density  Is  not 
a  complete  description  of  the  random  process.  In  a  sense 
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it  is  obtained  as  an  "averaging  process"  over  very  long 
lengths  of  terrain;  therefore,  a  "discrete  rut"  or  any 
kind  of  pothole  in  the  runway  is  "averaged  out".  In 
addition,  use  of  the  decaying  exponential  input  restricts 
analysis  to  a  linear  system  as  is  verified  in  Appendix  5; 
hence  both  the  (1-cos  v;t)  bump  and  the  decaying  exponential 
are  used  as  system  inputs  to  broaden  the  analysis. 
Illustrations  of  the  (1-cos  wt)  bump  and  runway  spectral 
density  complying  with  (K11-A-8862A)  design  criteria  for 
runway  roughness,  are  given  in  Figures  2  and  3. 


11 
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III.  The  OptlT.nl  Control  Approach 

Cost  Criterion 

An  optimization  criterion  for  an  active*"shock  and 
vibration  system  has  usually  been  to  minimize  the  absolute 
acceleration  of:,  the  sprung  mass,  while  holding  the  rela¬ 
tive  displacement  or  rattle  space  of  the  suspension  system 
below  some  predetermined  roaxium  value  (Ref  1,2,3  )• 
However,  large  aircraft  landing  gear  have  as  much  as  two 
to  three  feet  allowable  rattle  space  to  absorb  landing 
Impact  energy,  while  less  than  a  half  foot  of  displace¬ 
ment  Is  required  to- absorb  most  runway  Imposed  vibration 
energy,  the  displacement  requirement  being  directly 

-  •  »  •  "»  •  1.1.  .  v..  ...  .  _  ,1 . .  .  »  ..  ..  ..  -i  •  j. _ 

S'X  W  V  •-  l/V  l/‘  ’  VX  WV  |/X  VUUXVU  C*IU^X- A  V  • 

On  the  other  hand,  due  to  the  relatively  large  mass 
of  the  airplane,  there  Is  the  expected  requirement  of 
large  control  forces  from  the  active  element  or  actuator. 
It  is  therefore  felt  that  as  a  first  try,  the  optimization 
criterion  should  be  to  minimize  the  absolute  acceleration 
of  the  sprung  mass,  the  airplane,  while  at  the  same  time 
placing  a  constraint  on  the  control  force  instead  of 
rattle  .space.  In  other  words,  the  maximum  control  force  • 

f 

obtainable  is  the  constraint  resulting  from  physical 
limitations’ of  the  hydraulic  or  pneumatic  actuator  used. 

Vith  the  above  criterion  a  system  Is  visualized 
where  the  actuactor  would  function  only  in  a  rattle  space 

7  - 
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range  of  +  6  in.  allowing  passive  elements  to- set 

limitations  on  rattle  space  when  the  +  6  in.  ranre  is 

exceeded.-  The  active-passive  combination  system  would 

seem  especially  attractive  since  it-'  would  retain  the 
'  -  'i 

passive  system  in  case  of  active  system  failure. 

The  optimization  or  minimization  criteria  can  be 

expressed  as: 


rft  } 

Minimize  J"((i)  =  /  [  Z  -f  PU(iJ  J  e/f 

*  -m 


(4) 


/ 

/ 


where  U(t)  is  the  control  function  Fs(t)/Pl  and  zTs"  the 
absolute  acceleration  defined  by  equation  3. 

Actually,  there  are  a  number  of  mathematical 
formulations  that  could  be  used.  If  the  maximum  value  of 


Z  is  minimized  along  with  a  maximum  value  of  rattle  space 
XI,  the  solution  would  result  in  a  bang- bang  type  of 
controller  for  U(t)  (Ref  4,5).  If  the  criterion  is  to 
minimize  Z  in  minimum  time,  the  solution  would  again  lead 
to  a  bang- bang  type  of  controller  for  U(t).  However,  the 
bang- bang  controller  requires  comparably  larger  control 
forces  than  the  continuous  duty  controller,  and  it  often 
demonstrates  a  poor  ability  to  withstand  broad  band 

random  inputs;  thus,  the  bang-bang  controller  appears 

,  /  • 

undesireable  for  landing  gear  use  (Ref  3:14). 
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IiXcontrast,  the  integral  square  criterion  .of 
equation  (4)  offers  many  practical  advantages  for  the 
active  landing  gear  system.  The  resulting  system  obtained 
from  this  criterion  is  linear  and  therefore  easily  ’ 

analyzed,  and  a  constant  coefficient  feedback  system  can 
be  found  arid  can  usually  be  realized  by  an  active  system 
(Ref  3).  The  quadratic  criterion. is  also  advantageous  in 
that  a  global  minimum  is  assured.  In  other  words,  there 
is  only  one  control  function  U(t)  which  satisfies  the  above 
criterion  for  a  given  system  -with  a  specified  set  of 
parameters.  Finally,  the  integral  square  criteria  is 
required  when  the  equivalent  deterministic  input  is  used 
for  the  system  must  be  linear  cefore  mean  square  output 


VAl  UAA  fflfl  hp  f.VlPAT’O  t*  1  1  *1  'f  n>»£»r?  4  ♦.«  4  . 

'  **  a  —  - - -  - - -  ■*•*»•  w**  w  jl.  v/;u  fc*ii 

input  derived  from  runway  spectral  density  (Ref  Appendix  B). 
Of  course  the  criteria  in  equation  4  for  use  with  the 
equivalent  deterministic  input  can  be  stated: 


Minimize  T(u)  =  JE  f  Z  )  +  P  E  ( U?  (-tf)  (5) 

where  the  ?.  (  )  notation  denotes  "expected  value" . 

Theory 

Equations  4  and  5  establish  the  cost  criterion. 

To  continue  with  the  problem  formulation,  by  substituting 
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equation  2  and  3  In  equation  1,  the  system  equations  in 


state  space  notation  are 

€ 

Xi  =•  Xa.  (6) 

=  (l+R)JU(t)  -O3X1  -Q4X1 
—  &  7  (  !=(«-)  -  Xa)  (?) 

Xl  =  X  4-  (0) 


X't  —  R  Cj  U.  ( (r)  -f  Q  £,  Xi  -f  Q7  X* 

4-  r  (  tr/e'i  ^  Y-i  )  —  /oi 

‘  ' . .  '  J 

where  R,  g,  Q3,  $4,  4*5 ,  Q6,  Q7,  are  defined  in  the  list 
of  symbols. 

The  Hamiltonian  function  (Ref  10)  corresponding  to 
the  system  differential  equations  6-9  Is 

n 

H  -  Z  Pz  Xi 

x!-0 

(10) 

=  Po  (  x  2  +  P  u^it)  ) 

.  r?  V 


CGr/K-:/70~6 


-t  Px  [  (  f  +  A)  J  U  (O  -  Q3X1 
~  Ql  Xa  ~  QF  (  F^-Xb)] 


+  P3  XI 


t  P  4  [  ~  Rq  U  (i-)  +  Qi)i,  f  0  7  X  a 


+■  Qt  C  F  (+)  -  X  3  )  - 


vh^re  ’  x  1  r,  I'JiUort  the:  joint  variable.  It  corresponds 

to  the  L-’'ra:a'e  multiplier  In  the;  conventional  calculus 
of  variations  vhen  equality  constraints  are  present. 

The  system  differential  equations  of  the  adjoint 
variable  are  derived  as 


Px  =  - 


Ai± 

J  Xi 


X  -  1,3. 


and  the  trnnsversality  conditions  are  determined . by 

P/M  -  iJ=L  [X/f,)  ql  i  -  i) } 


ggc/se/70 -6 


where  ©  is  that  part  of  the  cost  function  which  is 
a  function  only  of  final  time  tf.  In  other  vfords,  in 
terms  of  the  gerieral  cost  formulation  of  the  Rolza 
problem,  the  cost  expressed  In  equation  4  contains  no 
function  ©  which  is  a  function  of  tf  only;  therefore, 
©  is  zero  (Ref  10:57). 

Finally,  for  a  minimization  problem,  };0  i  rs  a  const  a 
and  greater  than  zero  (Ref  10:64).  The  Hamiltonian  is 
homogeneous  in  P^*,  therefore  I© can  be  set  eqi:al  to  one. 
For  U(t)  to  be  optimal,  the  Hamiltonian  function  H  has 
to  be  a  minimum  with  respect  to  the  variable  U(t). 

A  necessary  condition  for  optimality  is  then 


d  H  =  O  U3) 

dU 

,  - 

•  -  .  • 

The  solution  of  equations  6-9  and  11-13  gives  the  optimal 
solution  for  U(t).  A  detailed  mathematical  formulation 
of  the  problem  can  be  found  in  Appendix  A. 


/ 
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IV .  Numerical  Techniques 


Conjugate  Grad l ent 

The  optimal  control  numerical  technique  used  to 
solve  the  problem  of  this  study  Is  the  conjugate  gradient 
technique  of  Lasdon,  Hitter  and  Warren  (Ref  8).  The 
technique  is  directly  applicable  only  to  unconstrained 
problems ,  but  the  constrained  problem  can  be  converted  to 
an  unconstrained  problem  through  the  use  of  a  penalty 
function  as  Is  done  In  equations  4  and  5*  The  uncon¬ 
strained  problem  under  consideration  Is  to  minimize  the 

absolute  acceleration  Z.  At  the  same  time,  a  penalty 

•  # 

must  be  paid  In  acceleration  Z  to  constrain  the  control 
force  U(t)  bcloi’?  some  predetermined  value  set  by  the 
physical  limitations  of  the  actuator. 

The  conjugate  gradient  method  Is  a  direct  method. 

In  other  words,  for  an  optimal  solution,  a  search  is 
made  In  the  direction  which  directly  minimizes  the  cost 
function.  The  gradient  trajectory  (Ref  equation  13) ,  its 
norm  and  the  actual  search  direction  are  the' only  values 
which  require  computer  storage.  The  search  directions 
are  formed  from  past  and  present  values  of  the  cost  and 
Its  gradient.  Then  successive  points  are  determined  by 
linear  minimization  along  the  search  directions  which  are 
always  directions  of  descent.  With  search  directions 
always  descending,  the -con Jugate  gradient  tends  to  converge 
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even  from  poor  approximations  to  the  minimum. 

Algorithm 

The  conjugate  gradient  algorithm  requires  the 
computation  of  the  gradient  trajectory.  Letting  UQ(t) 
be  the  first  of  approximations  to  the  optimal  control 
Ui (t ),  1  =  1,2  etc.,  then  the  corresponding  gradient 
G0(u)  Is  computed  by  solving  the  state  equations  6-9 
forward  in  time,,  the  adjoint  equations  11  backward  In  time 
and  obtaining  the  gradient  dii/dU  from  equation  13. 

The  algorithm  proceeds  as  follows : 

1.  Guess  U0(t),  an  arbitrary  selection  of  UD(t)  =  o 
was  made 

<■*»  T.  1  .  —  i  -  .  «  -  I  .  ...  **  ....1  -  ."I  __T»  I  .  .  i 

*-  •  O.  I  k  tAW  H'-4  VlOUkl  1  OA  IIVUU  l*U  joauk 

equations  backwards 

3.  Compute  gradient  0o  =  G0(uo)  from  equation  13 
Compute  search  direction  SQ  =  _g0 

5.  Choose  oc  =  oc  to  minimize  J(Ui+ cC^S^ ) 

6.  Ui+i  = 

7.  G1+1  =  G(U^+^) 

8.  =  /  (Gl+1  »^i+l  )/(Gl  *G1 ) 

9.  Si+1  =  -G1+1  +  l^S, 

where  (G^.G^)  = 

10.  Set  1  =  1  +  1,  go  to  5 


2.0 
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Steps  5-10  are  repeated  Integrating  the  state  equations 
forward  and  adjoint  equations  backward  for  each  com¬ 
putation  of  the  gradient,  until  the  change.  In  cost  is 
ne-llgible. 

Alpha  Search 

The  alpha  value  indicated  in  step  5  of  the  algorithm 
is  determined  by  initially  assigning  a  small  value  -  about 
.1  -  to  alpha  and  then  checking  the  inner  product  Gj[+i  3^. 
If  the  sign  of  the  inner  product  is  positive  with  the 
Initial  alpha  value,  a  tenth  of  the  initial  value  is 
repeatedly  token  until  a  value  of  alpha  is  found 
which  gives  s  negative  inner  product.  Then  the  alpha 
value  is  repeatedly  doubled  until  the  inner  product  is 

r*  ^  4*  1  **A  f***  V>  4  A  4  v%  ♦-  1  r.  4*  *  1  4*t»  #%•«»  i«M«k4  ^  4.  —  f  • 

*  -  -  -  —  -  —  -  -  -  •  -  r  -  -  •  •  -  *»•  '  wvb  W  wfcO  VV»4.  014.110 

the  value  of  alpha  which  gives  a  zero  inner  product.  The 

.  •  • 

•  ^  * 

inner  product  or  directional  derivative  is  the  slope., of 

•  :  •  V  ...  \ 

the  cost  with  respect  to  alpha;  therefore,  the  alpha 
value  with  a  corresponding  zero  inner  product  is  the 
alpha  value  which  gives  the  lowest  cost  in  that  search 
direction. 

Some  alpha  searches  use  one  over  the  norm  of  the 
search  direction  3  as  a  guide  for  the  first  initial  alpha 
guess  (Ref  o).  However  after  trying  this,  it  was  found 
from  running  a  number  of  sample  problems  (Eef  8:136)  and 
the  dynamics  of  the  thesis  problem  that  a  guess  in  the 
range  of  0.001  to  0.1  nearly  always  gave  a  negative 


\ 
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inner  product  and  was  prenerally  better  than  that  of  one 
over  the  norm  of  the  search  direction.  A  better  alpha.' 
pruess,  of  course  means  a  savings  in  computer  iteration 
time.  The  computer  program  can  be  found  in  -\pperv:  2. 


ggc/se/70-6 

V.  Results 

All- results  are  based  on  the  following  arbitrarily 

selected  parameters  for  the  dynamic  model. 

I  Cs  100  lb  sec/ft 

\  PI  100,000  lb 

\  . 
v- — 1,000  lb 

„  -  Ks  200,000  lb/ft  - 

Kt  1200,000  lb/ft  . 

The  response  of  the  dynamic  model  to  the  (1  cos  v/t.) 
bump  input  is  illustrated  in  Figures  4-17.  The  plots  shov: 

9  9 

the  control  force  U(t)  required  and  the  acceleration  Z 
piroduced  when  the  weighting  factor  F  is  set  equal  to 

0,  1,  10  and  infinity.  The  value  of  P  equal  to  zero 

•  « 

illustrates  the  virtually  zero  acceleration  of  the  aircraft 
produced  when  there  is  no  limit  placed  on  the  control  force 
--  iKt).  With  F  equal  to  one,  equal  weighting  is  then  given 
in  the  cost  function  to  minimize  both  the  control  force 
U(t)  and  the  acceleration  2  (see  Appendix  A).  With  a  P 
value  of  10,  control  is  weighted  even  more  until  finally  . 
with  F  equal  to  infinity,  the  control  force  is  then  zero' 
and  the  resultant  acceleration  Is  the  acceleration 

v 

produced  by  a  completely  passive  system. 

Figures  4-1?  cover  £he  frequency  range  of  .5  to  JZ  Hz 
geometrically  by  doubling  the  frequency  of  the  argument  of 
the  cosine  function  beginning  at  .5  Hz.  The  bump  height 
or  amplitude  of  the  (1-cos  wt)  function  is  2.5  in. 

1  ■'  ‘ 
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corresponding  to  the  paved  runways  in  Figure  2.  Th'e 
figures  are  arranged  in  pairs  such  that  the  control  force 
required  for  the  four  values  of  P  for  a  particular 
frequency  bump  is  given,  and  then  on  the  following  page, 
the  corresponding  acceleration  is  given  for  the  same  values 
of  F  and  same  bump  frequency . 

One  observation  that  can  be  made  is  the  increase  in 

control  force  required  as  the  bump  frequency  increases. 

»• 

This  is  reasonable,  since  the  faster  the  aircraft  moves 
i 


over  the  bump,  'the  faster  the  actuator  must  contract 

I 

the  landing  gear  and  then  extend  the  landing  gear  to 
maintain  the  FI  force  at  all  times.  The  higher  aircraft 
velocity  also  accounts  for  the  increase  in  force  required 


t«t4  4-la  /-*■*»***«->  r-  r\  A  v»  *>-. 


,  ......  „  — _ J  - ±.-\. 

»  »  WUItU  l/i 


■T*  J - 1  -.1  • 

X  XUOli. 


of  the  bump  'because  again  as  the  aircraft  moves  over 

the  bump,  the  aptuator  must  extend  the  landing  gear  faster 

to  maintain  the  '^aircraft  at  a  fixed  level. 

A  second  observation  is  the  waviness  in  the  plots 

at  the  lower  frequencies  with  the  weighting  factor  P 

equal  to  one.  The  weighting  of  P  equal  to  one  tends  to 

make  the  cost  function  work  against  Itself  because  it-  is 
1 

•  * 

being  asked  to  find  the  U(t)  to  minimize  Z  but  at  the  same 
time  keep  U(t)  small.  These  two  values  of  U(t)  of  course 
a r  ^  opposite  ends  of  the  scales.  In  addition  and  even 
more  important ,  the  dominant  resonant  frequency  of  the 
system  is  located  at  1.2  Ilz  which  docs  not  permit  the 
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computer  program  to  converge  to  ns  low  a  cost  near  1  ;.z 

/  . 

as  reached  at  other  bump  frequencies; 

A  third  observation  of  course  is  the  marked  improve¬ 
ment  of  the  active  system  over  the  passive  system.  sith 
P  equal  to  infinity,  there  is  considerable  vertical 
acceleration  to  the  aircraft,  but  v:  1th  F  equal  to  zero,  . 
the  acceleration  is  virtually  zero  meaning  that  «s  long 
as  the  actuator  can  produce  the  force  required  there  will 
be  little-or  no  acceleration  to  the  aircraft. 

For  the  decaying  exponential  input,  a  plot  of  the 
root  mean  square  (RMS)  values  of  airplane  acceleration 
and  actuator  control  force  against  the  weighting  factor  ;  . 
is  given  in  Figure  18.  The  RH3  values  are  obtained  by 
takluv  f.Vi r»  snuara  rnnt;  n,!*  t.hA  intAfi's  I  ufmaro  ^rtO»A  n+- 
acceleration' Z  and  the  square  root  of  the  Integral  square 
value  of  control  force  U(t).  Ha thematically  it  means 
evaluating  equation  4  as  two  separate  integrals  without 
the  weighting  factor  P. 

The*  RMS  values  for  the  acceleration  and  control  for 
F  equal  to  infinity  are  not  shown  on  the  plot  as  in  the 
(1-cos  wt)  plots  because  the  value  of  acceleration  for  the 
totally  passive  system  would  be  off  the  paper.  The  RMS 
value  of  acceleration  at  P  equal  to  infinity  is  .54-3 g 
and  the  HH3  control  for  this  value  of  P  Is  .000169  x  PI  or 
virtually  zero. 

For  the  RFS  plot,  the  decaying  exponjetial  input  is 

39  " 


WEIGHTING  Fr.CTGB 


GGC/S E/70-6 


derived  from  the  runway  spectral  density  of  the  form 
(2ft  \//Ac)x)  .  The  A  is  a  parameter  equal 

to  10“5  for  most  runways  (Ref  1:49),  and  V  is  the  velocity 
of  the  aircraft  which  is  set  equal  to  200  ft/sec  to 
represent  a  maximum  take  off  velocity.  The  JTTV/Ao 
term  represents  the  break  frequency  at  vrhich  the  runway 
spectra  level,  off.  The  runvray  spectra,  of  course  must 
level:  off  due  to  the  finite  height  of  runways.  Using  an 
?.VS  value  of  runvray  height  of  .5  ft  (Ref  1:14),  and  the 
above  spectral  density  form,  the  value  of  the  wave  length 
•^0  is  found  to  he  5*10*  ft  from  equation  33  in  Appendix  B. 
The  above  values  of  A, 7,  and  A®  are  then  substituted  In 
equation  36  of  Appendix  B  to  obtain  the  exponential  input 

\Zii1T  A  V  EXP  ,3  IT  \/  ~X  /  Ae^  need  fnv  4-V»o  i?»  nin»o  1  ft 

'  •  •  -  —  - 

plot . 

The  plot  reveals  the  trade-off  between  RMS  accelera- 

09  J 

tlon  Z  and  control  force  U(t).  The  maximum  value  of  RES 
control  force  required  is  of  course  where  P  equals  zero. 
This  value  is  approximately  7200  lb.  In  other  words  the 
7200  lb  is  the  standard  deviation  of  U(t)  or  c fu. 

Therefore  given  a  normal  distribution,  3°^  oi*  a  21600  lb 
force  would  be  the  maximum  force  required  from  the  actuator 

»  a 

to  ensure  that  the  probability  of  actuator  ‘saturation 
•would  be  less  then  .3  percent.  Like  the  (1-cos  wt)  input 
plots,  the  RMS  plot  also  reveals  the  better  vibration 

t  t  * 

Isolation  performance  of  the  active  system  over  the 
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passive /system.  As  the  weighting  factor  P  increases,  t^e 
acceleration  becomes  as  high  as  .543s  with  the  completely 
passive  system  at  P  equal  to  infinity,  but  with  the  active 
controller  in  the  system  the  acceleration  can  be  kept  beloi 
an  RMS  value  of  .07g.  ^ 

Numerical  Aspects  * 

For  the  (1-cos  wt)  input  the  computer  program  con-, 
verged  in  4  to  5  iterations  in  approximately  3  minutes 
time.  One  Exception  was  the  convergence  rate  around  the 
system  resonant  frequency  of  1.2  Hz.  Near  this  frequency, 
a  U(.t)  guess  of  zero  or  other  constant  value  guess  would 
require  15  to  20  minutes  computer  time  to  converge  to 
approximately  the  same  cost  values  obtained  at  frequencies 

nwnv  frnw  tho  Con 

vergence  near  the  resonant  frequency  was  as  fast  as  that 
of  other  frequencies  if  the  solution  for  U(t)  found  away 
from  the  resonant  frequency  was  used  for  frequencies  near 
the  resonant  frequency. 

With  this  procedure,  the  number  of  iterations  arid 
the  time  to  converge  were  reduced  to  5  iterations  and 
3  minutes  computer  time  respectively,  but  the  cost  values 
obtained  near  the  resonant  frequency  were  still  not  quite 
as  low  as  those  of  other  frequencies,  as  Is  evident-  by  the 
waviness  of  the  1,2,  and  4  Hz  plots  of  the  (1-cos  wt)  bump 
input.  ~At  these  frequency  Inputs  the  accuracy  of  the 
values  of  U(t)  and  Z  fell  from  3  significant  figures  to 
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about  one  significant  figure  or  a  "ball  park"  value. 
Possibly  In  the  resonant  frequency  area  a  "second 
variational"  numerical  technique  would  give  better 
solutions  for  U(t)  and  Z  using  the  gradient  solution  for 
U(t)  as  an  estimate  (Kef  10). 

Another  significant  factor  which  effected  the 

convergence  rate  was  the  definition  of  U(t).  At  first 

U(t)  was  defined  F(t)/?1.  This  resulted  in  veirh tings 
♦  ■ 

in  the  dynamics  of  the  ^odel  with  ratios  ' of  1  to  12g ,  that 
Is,  the  I'relghting  of  U(t)  would  be  1  while  that  of  another 
term  might  be  12g.  In  contrast,  by  defining  U ( t )  to  be 
F(t)/Fl  (see  Appendix  A)  the  highest  weighting  ratio 
between  any  two  terms  became  ig  to  12g."  The  end  result 
of  this  definition  was  that  the  number  of  iterations  and 
the  time  for  convergence  was  cut  in  half.  For  example, 
for  the  (1-cos  v;t)  inputs  the  average  number  of  iterations 

for  convergence  was  reduced  from  10  to  5. 

»  * 

The  convergence  time  for  the  decaying  exponential 
Input  was  considerably  longer  than  that  of  the  (1-cos  wt) 
input.  It  generally  fell  in  the  range  of  20  to  30  minutes, 
for  various  values  of  the  weighting  factor  P.  The  longer 
convergence  time  was  expected  since  the  decaying  expon¬ 
ential  input  had  a  very  slow  decay  rate'  which  required  .the 
mod^—dyTurmTcs” to  -be  solved  over  a  comparably  greater  real 
time  than  that  required  by, the  (1-cos  wt)  input.  However, 
It  was  still  felt  that  maybe  the  4th  order  Runge-Xutta 
integration  procedure  In  Appendix  C  was  not 

c 
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sufficiently  accurate,  there-by  causing  a  longer  con- 

i. 

vergence  time.  Consequently,  another  gradient  prorr.xm  was 
written  using  a  prewritten  Integration  routine  and  a 
function  minimization  routine  already  on  hand .  The 
function  minimization  routine  was  also  modeled  after  the 
conjugate  gradient  method  as  the  one  in  -vopcndlx  C. 

The  second  program  required  more  time  per  iteration, 

j 

while  the  cost  decrease  was  greater  per  iteration.  The 

•  I 

Appendix  &  program  required  less  time  per  iteration,  but 

the  decrease  per  Iteration  was  not  as  great.  However 

both  programs  converged  to  almost,  identical  costs  given 

equal  time.  For  example,  with  the  weighting  ?  set  equal 

'to  .1  in  each  program,  the  program  in  Appendix  C  converged 

to  a  cost  value  of  ,470  x  10-2  while  the  second  n-rog-rnrn 

converged  to  a  cost  value  of  ,468  times  10”2. 

Convergence  Criteria 

The  value  of  the  cost,  the  gradient  and  the  Kamil-  * 

m  * 

tonian  function  can  De  used  to  check  for  convergence  of  the 
numerical  problem.  For  the  results  of  this  study,  the 
values  of  cost  and  gradient  were  checked  as  well  as  the 

values  of  X2  and  X4.  In  checking  the  dynamics  and  the 

\ 

model  of  Figure ;1,  It  can  be  geen  that  if  the  weighting 
factor  P  is  zero,'  the  value  of  the  relative  velocity  X2 
between  Ml  and  K2  is  equal  and  opposite  in  sign  to  X4,  the 
velocity  of  M2  relative  to  the  fixed  reference.  This  must 
be  so  for  the  airplane  Ml  to  ride  perfectly  level  over  a 
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/ 

/ 

/ 

bump  or  depression.  For  the  weighting  factor  P  equal  to 
zero,  these  values  agreed  in  magnitude  to  within  three 
and  four  significant  figures.  Finally  the  gradient  values 

p  °  o 

for  convergence  were  In  the.  ranjre  of  10  for  the  discrete 

/  *  / 

time  points  and  the  chanre  In  cost  value  between  Iterations 
for  convergence  was  set  at  less  than  .5  percent.  If  the 
computer  program  was  left  to  Iterate  after  the  above 
criteria  had  been  net,  no  significant  improvement  In 
cost  was  gained.  / 

v..,. 

/ 

/ 
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VI.  Conclusion 

The  response  of  the  active  landing  gear  system  to 
both  the  (1-cos  v:t)  Input  and  the  decaying  exponential 
Input  shows  a  marked  improvement  in  vibration  isolation 

performance  over  the  passive  system.  \ To  achieve  the 

l 

active  system,  'the  actuator'  is  selected  and  its  maximum 
values  of  force  output  and  response  rate  arc  noted. 

'With  these  values,  the  weighting  factor  F  in  the  cost 
criterion  is  chosen  such  that  the  system  response  requires 
less  thqn  or  the 'equivalent  of  the  maximum  force  and 
response  limitations  of  the  actuator.  The  designer  can 

then  weigh  the  vibration  isolation  performance  to  be 

* 

a  ■  _ 

gained  by  the  active  landing  gear  system; 

■  ,  ■» 

The  study  provides  "open  loop"  solutions  for  U(t) 
for  both  kinds  of  inputs;  however,  the  solution  for  U(t) 
with  the  decaying  exponential  has .thb  advantage  since  *  • 

.  -  4  <  ’  V. 

the  cost  function  is  the  integral  square  criterion  with 

.  V 

the  final  time  tf  approaching  infinity.  -With  this'  criteria 
a  constant  coefficient  feedback  controller  can  be  obtained 
from  the  open  loop  controller  U(t)  (Ref  10, 9i). 
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VII.  Recommendations 


/ 


During  the  course  of  this  study,  a  number  of) 

\ 

_  V  ; 

Investigations  for  ^further  study  have  been  nur nested . 
The  following  Is  a  list  pf  possible  experiments ; that 
might  prove  valuable  for  the  active  control  landing 'pen: 

1.  The  rattle  space  could  be  constrained  and 

\ 

*  ..  t 

made  part  of  the  cost  function  along  vfith 
the  control  fore*-.  Uj(t)..  " 

2.  Nonlinear  passive • elements  could  be  used 
to  more  closely  simulate  present . landing 


gear  suspensions. 

3.  Improve  the  dynamic  model  of  this  study  by' 

providing  for  wheel  diameter  size  effects  and 
permitting  wheel  hop  for  very  rough  landing 
surfaces .  * 

4.  Investigate  model  response  with  various  values 
of  Ks  and  Cs. 
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with  tf  equal  to  the  period  of  the  cosine  function  for  the 
(1-  cos  v:t)  input  and  equal  to  infinity  for  the  decaying 
exponential  input. 

Kow,  in  order  to  evaluate  equation  1 6  in  the  same 
integration  scheme  as  equations  6-9,  a  new  state  variable 

h<) 


n-c/^/7  o.f. 


Is  Introduced  such  that 


Xs  =  l  (utt)  -Oi  X,  -QaXx-i) 

+  p  uact)  J 


2 


vrhere 


Xr(o)  ~  c 


Then  equation  16  becoir.es 


X(u)  -  X5  (tf) 


and  the  system  differential  equations  are 


Xi  —  Xa, 


-  Q?  (FCt)  -X3  ) 


Xl 

X+ 


=  X« 


—  R  0  L l(t)  +  Qi  X i  +  Q7  Xa 

-t-  Qs  (F  (i)  -  X.  3  )  ~  Q 


(i°.; 


Xa  =  (H-R)flU(t)  -Q3Xi-04Xa, 


(19) 


X$  —  Cu(t)  —  0 1  Xi  ”"Qa.  Xa  "I )  -f*  P  U.  (t) 


50 


'  ’/  Y'O— G 


-51 


P,  =  Q?  Pa  -  Q&  P4 

+  2  Qi  (uCt)  -  Oi  X*  ~~  fa  —  I  ) 

Pa  =  -  P.  +  Q4  Pj.  “  Qt  P 4 

t  2  Qa  CufO  ~  Qi  /i  —  0=  X  =  I  ) 

P 3  —  Q 5"  (  P'r  ~  Pa  ^  (22) 

P4  =  -Pi  ■ 

Pi'  =  o 

# 

id.  frorp  eqvuH  on  12,  He  o relieve r 3 c.3. it y  condii-lcms 

*G 

P.Ui)  =  o 
Pa.  (tf)  -  O 

pj  ( tyV  —  0  (25) 

P4Ct>)  =  O 

*■* 

PH)  =  i 


52 


?s 


T  l  T,-»  I  *  '  t*~*  ■**' 


‘s  -vw  v'..^ 


p?  (tv)  of  c.ai^ti^HS  28  r:-ii  29  resi  eov 

•o'. \ ; -oci:io:t  to  ive  P  5(t)  =  I- 
"Silent  of  cr-vcit 13  becones 


dH  =  C 1 1 r ) q  Pi  -  r q  P+ 
du 

t  2  (  UfO  -  Qi  Xl  -Qa  ) 
t  x  P  U  (-t)  (2-:) 
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Appendix  B 

Derivation  of  Input  F(t )  From 
Spectral  Density 


The  effects  of  runway  uneveness  can  be  represented 
by  a  stationary  broad  band  random  process  (Ref  7). 
Referring  to  Figure  3,  runway  height  Is  commonly  described 
in  terms  of  power  spectral  density.  The  spaclal  frequency 
_n_  (rad/ft)  can  be  transformed  to  the  time  frequency  LJ 
(rad/sec)  by  noting  that  the  distance  based  spectrum 

Is  related  to  the  time  based  spectrum  <^xx  (iv) 


<Mw>  = 


-O-  =  m 

V 


(26) 


and  V  is  the  velocity  of  the  airplane.  From  Figure 3 
the  power  spectrum  can  be  approximated  by  the  straight 
line  fit 


$xx(-n)  =  .A. 

.nt 

or  using  equation  25  and  26 

<kxM 


(27) 


(28) 
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A  Is  a  pawretfir  which  is  equal  to  10”^  for  most  run¬ 
ways  (Ref  1  :k?).  However  equation  27  becomes  infinite  rt 
the  low  frequency  end  where -as  elevation  spectra  must 
level  off  due  to  the  finite  height  of  runways.  This  can 
be  C1.  npensated  for  by  modifying  to  the  form; 


^XX  “  - =-^- 


-TC  -h/J£L\ 


VAo  / 

and  again  converting'  to  the  time  based  spectrum 

4>xx(w)  =  av 


(JO, 


\  A< 

where  A®  was  chosen  to  be  5 '  X  10^  ft.  (Ref  1:47). 

The  nhnpo  ponnf  \  on  wf»»ouonfr  m  O  <««e  1 

*  -  —  **  r-  -  -  -  —  -*•  — 

density  or  the  time  based  spectrum  for  an  e rootle’ random 
process.  Therefore  4>*xM  is  the  fourier  transform 
of  the  autocorrelation  Rxx<rj 


.  f  *° 

Rxx(r')  =  4.  /  <f>xx  (s)  es*ds 

+  L\  <- 


OD 


Ss/» 


The  objective  is  to  find  an  F(t)  so  that  one  can 
minimize  the  nean  square  value  in  equation  5.  The  mean 
square  value  pf  F(t)  is  Rxx(o) 

or 
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Fa  =  Rxxfo) 


_L  /  (pxx(S)  <Js 

> 


Using  Parseval 1  s  the  or? ,  the  squ'-r*-  veV.)--  '-an 

be  vrritten 

F(s)F(-S)dS 

oo 


where  F(s)  is  the  four! er  transform  of  ?(t)  ( ?.e f  6). 

Kow  if  (j)  X  xfs)  is  Ku.lt i pi i O'"  by  JlTT/4iTT  and  the 
denominator  term  is  taXen  outside  the  integral  in  equal- io: 
32,  then  comparing  equations  32  and  33 


F(S)F(-S)  =  <£xxfs)  =' 


art  av 


(3*0 


The  time  function  F(t)  is  equal  to  the  inverse  transform 
of  F(s),  or 


F(t)  =  F'1 


/atTAV 

s  +  3  ttv 

Ac 


F  (t) 


=  u«3  van  av  e 


-a  rr  v-t 


(35) 


(36) 
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decaying;  exponential  Input  as  a  measure  of  rrand otsc-ss 

of. runway .height,  and  by  ensuring  a  linear  system,  root 

mean  square  system  output  values- of  force,  displacement, 

§ 

velocity  et <?•«'•  can  be  obtained  by  computing’  the-  square 

*  x/ 

root  of  the  integral  square  value  of  the  system  output 

r.  r 

force;  displacement:  and  velocity  respectively. 

**  % 

«  ' 
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.  Appendix  C 

Digital  Computer  Program 

The  IB”  7090  digital  computer  was  used,  to  make  the 
calculations  found  in  previous  sections  of  this  study. 

The  program  as  Illustrated  in  the  following  pages,  is 
arrahged  for  the  (1-cos  wt)  input,  ilowever  by  not  storing 
all  the  states  and  F(t),  more  of  the  necessary  computer 
storage  space  can  be  gained  for  use  with  the  decaying 
exponential  Input.  The  list  of  symbols  in  the  prefatory 
of  this  study  also  applies  to  the  computer  program. 
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$  I B JOB 
SIBFTC  MAIN 

COMMON/GOOD/X(5*.?OOl  )  »XX<5«5>*DXC5«5>  *DH<5)  * 

1G (2*2001 ) *S(2001 ) *  U  ( 2  * 200 1 ) *F(2001 ) « KCOUN* GRAV 
2C0ST «NNZP«LB0  *M«  P « 01  *  02  *03  *04  «Q5 « 06«C7«R 
DIMENSION  ACC (501) 

C  H  15  the  DISCRETE  TIVF  INTERVAL 
C  IT  IS  THE  NUMBER  OF  ITERATIONS 
C  KCOUN  IS  THE  NUMBER  OF  DISCRETE  TIME  POINTS 
FREQ  =  8*0 
KCOUN  =  401 
..RCflNM  =  KCOUN  -  1 
H  =  1 «0/(FREO*RCOUN) 

KSTER  =10 

C  GUESS  FI PST  CONTROL 
DO  16  I  s  1 i KCOUN 
16  U( 1 « I  1  =  0.0 

IS  IT  ■  0  ’  /  •  ‘  f 

TEM  =  1.FC5 
NNZP  =1  ~ 

NN  =  -1 
LBO  =  -1 
BETA  e  0.0 
PI  12  =  6*26318 
GRAV  =  32.1725 
F i  =  innnnn. 

F 2  *  1 OOO* 

XM1  =r FI /GRAV 
XM2  =  F2/GRAV 

XCS  =  100.0 
XKS  =  200000* 

XKT  =  1200000* 

P  =  0.1 
01  =  XKS/F1 
02  =  XCS/F1 

03  =  XKS*GRAV*(F1+F2)/'(F1*F2>  ^ 

04  =  XCS*GRAV*(F1+F2)/(F1*F2> 

05  =  XKT  *GRA V/F2 
06  =  XKS*GRAV/F2 
07  =  *CS*GRAV/F2 
R  =  X&I1/XM2 

FS  *  ( SORT ( XK  S/XM 1  ) ) /P I  I  2 
FT  =  (  SORT ( XKT /XM2 ) ) /P I  I  2 
WRITE  (6*801)  FS»FT 
B01  FORMAT  ( 2X« 2TE1 5* 8 « 4X > ) 

TIM  =  0.0 

DO  1  I  =  1. KCOUN 

SCI )  =  0.0 

F ( I )  =  <1.0  -  COS( PI  I 2*FREQ*TIM) )*0. 10416667 
TIM  =  TIM  +  H 
DO  1  J  =  1*5 
1  XCJ*I)  =  0.0 


60 


ggc/27./70-6 


C  INITIAL  CONDITIONS 


X  (  3 . 1 

) 

=  -(F1+F2)/XKT 

X(  1  »  1 

>  =  -Fl/XKS 

M  =  1 

XX  (  1  t 

1  ) 

=  X(1  .1  ) 

XX  (2. 

1  ) 

=  X  ( 2  ♦  1  ) 

XX  (3. 

1  ) 

=  X  (  5  «  1  ) 

XX  (  4  « 

1  ) 

=  X  (  4  4  1  ) 

XX  (  5  * 

1  ) 

C  X  (5  «  1  ) 

C  INTEGRATE  STATE  EQUATIONS  FORWARD 
CALL  DFR1  (H.  1  ) 

DO  2  !  =  2 .KCOUN 
M  =  I 

CALL  PUNGE  (H«l«5) 

DO  2  J  =  1*5 

x(j«n  =  xx(j*2) 

DX(J«  1  )  =  DX( J«2) 

2  XX ( J . 1 )  =  XX  <  J  * 2  ) 

C  INTEGRATE  COST  ATE  EQUATIONS  BACKWARDS 
‘  V.  =  KCOUN 
H  s  -H 

C  FINAL  CONDITIONS  ON  COSTATFS 
XX(1«1)  =  0.0 

xx(2.n  =  o.o 

XX (3*1)  =  0.0 

VV  /  —  «->  .  n 

CALL  DER2  (H«l) 

DO  3  I  =  2 . KCOUN 
CALL  RUNGE  (H*l,4) 

DO  3  J  =  1 « 4 
DX< J. 1 )  =  DX( J.2) 

3  XX( J* 1  )  =  XXI J, 2) 

H  =  -H 

C  compute  COST 

COST  =  X  (  5  ♦  KCC-UN ) 

IF  (  NNZO  )  5.5.6 

C  ALPHA  SFARCH  AND  NEW  CONTPOL 

5  CALL  ALPHA  (H*KOK) 

GO  TO  4 

6  CONTINUE 

DO  26  I  =  1 .KCOUN 

26  ACC  (I)  =  U(1*J)  -  01*X(1.1)  -  02*X  (2.1)  -  1.0 
WRITE  (6.200)  IT.FREO 

200  FORMAT  (//,  5X «  12H ITERATION  =  «  16*1  OX* 

1 1 2HFREOUCNCY  =  .E1E.P.//) 

IT  a  IT  +  1 
WRITE  (6*201 ) 

201  F  CRM  AT  (  5X  »  4  HT  I  f-*E  .  ]  py ,  RHG°AD  I  ENT  *  1  1X.9HX1  (STATE)  « 

1 ox *PKX3( STATE) 

2»?X« 1 1HU1 (CONTROL) . EX. ?HX2( STATE) «  OX . 9HX4 ( STATE ) ) 
TIM  =  0.0 
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TSTEP  =  KSTEP 

DO  20?.  I  =  1 .KCOUN  4  KST-P 

VPITE<6.20  2>  TIM  .  G( 1 . 1) .X( 1 4  I ) .X<3. I >.U< 1 . I) . 
1X<2. 1 ) 4 X ( 4  *  I ) 

202  F OPM AT  (7(F1  5.B.4X  )  1 

2 03  TIM  =  TIM  +  H*TSTEP 
'  WRITE  ( 6  4  503  ) 

503  FORMAT  </////) 

WRITE  (6.204)  COST 

204  FORMAT (//.5X.6HCOST  =  4  E15.8  ) 

C  COMPUTE  SUM  OF  ABSOUiT'f 

C  VALUE  OF  GRADIENT  POINTS 
DIFF  =  C. 

DO  7  I  =  J  .KCOUN 
DFV  =  ADS  <G< 1 1 II) 

7  DIFF  =  DIFF+  DFV 
V.'PITF  (64205)  DIFF 

205  FORMAT  (14HSUM  OF  DIFF.  =4515.04//) 

C  CHECK  COST  AND  SUM  OF  ABSOLUTE 

C  VALUE  OF  GRADIENT  POINTS 

IF  (IT  .LT.  5)  GO  TO  706 

IF  (  DIFF  -  0.5  )  7014  701  4  702 

702  IF  (ABS( ITEM-COST) /COST)  -  .005)  701,7014706 
706  TFM  =  COST 

I F  ( NN )  1 0  4 1 O  4  3 

C  v. 1  t't_  .  M 

8  PP02  =  0.0 
PR04  =0.0 

DO  9  1  =  1  » KCOUN 

PROl  =  G ( 1 4 1 )  **2 
PR02  =  PP02  +  PROl 
PR03  =  G ( 2.  I )**2 

9  PR04  =  PPC4  +  PP03 

BFTA  =  PROP/PR04 

C  STORE  CONTROL  AND  GRADIENT 
C  COMPUTE  NEW  ^SEARCH  DIRECTION 

10  DO  11  I  =  1 4 KCOUN 
U ( 2 4  I  )  =  U(1 4 I) 

G <2f I)  =  G( 1  »  1  ) 

11  S(I)  =  — G (241)  +  BETA*S(I) 

NN  =  1 

NNZP  =  -1 
GO  TO  4 

701  WRITE  (6.206) 

206  FORMAT ( 1  OH#***#***** 4  21HLAST  TRAJ.  IS  OPTIMAL 
1  10H*********-*)  . 

WRITE  (7.311)  (  U (  1 4  I  )  4  I  =  1 4 KCOUN ) 

WRITE  (7.311)  (  ACC ( I  1 4  I  =  l.KCOUN) 

811  FORMAT  (5F15.S) 

IF  (FPF0-5.0)  31.31.32 
31  IF  ^tFPEO-1 .5)  33.33434 
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33  IF  ( FRFO— 0.6 )  37  «  37 «  36 
3?  re-O  =  FPEO  /2.0 

GO  TO  35 

34  FPEQ  =  FREQ  -  1.0 
GO  TO  35 

34  FREO  =  FREG/2.0 

35  H  =  1  •0/’<FREO*RCOUN) 

GO  TO  15 

37  CONTINUE 
STOP 

END  - ~ 

SlBFTC  SUBA 

SUBROUTINE  ALPHA  (HiKOK) 

C0MM0N/GOOD/X(5«2001 ) *  XX < 5 *5 ) «DX C 54 5 >  «DH<5>  « 
lG(2«2O01  1  4 S (2001  )  » 1X24  2001  ) 4r (2001  >  «KCOUNiGPAV* 
2CCST«NN2P«LB0*M*P  «C1 4 02 t 03 *  04 , QE « OS  *  07* R 
DIMENSION  A ( 60 ) «C0SS(60>  *PR04(60) 

30  IF  (  LEO)  10«  10 4  20 
C  ALPHA  COMPUTATION 
10  J  =  1 

AAA  =  2.0 

C  MAKE  FI PST  GUESS  ON  ALPHA 
AA  =  0.01 
A(J)  =  AA 

C  UPDATE  CONTROL 4  RETURN  AND  INTEGRATE 
DO  1  T  =  p.k'rnitM 

1  U(1«I)  =  U  (  2  4  I  )  +  A  (  J }  *S  (  1 1 
NNN  =  0 

LRO  =  1 
LL  =  -1 
RETURN 

20  PR04CJ)  =  0.0 
COSS(J)  =  COST 
NNN  =  NNN  +1 

IF  (NNN-EO  )  23424124 
C  COMPUTE  DIRECTIONAL  DERIVATIVE 
C  OR  I NT ERPPODUCT  G*S 
23  DO  21=  l«KCOUN 

PPOl  =  G< 1 4  I )  *  S< I) 

2  PR04IJ)  =  PR04IJ)  +  PRO  1 

V;RITF  (6*301)  PR04  ( J)  «A<  J)  *COSS<  J) 

801  FORMAT  ( 2X  < 1 C  HPR04 (J)  =  » E 1 5 . 8 1 4X *7HA ( J )  =  *E15.8. 
1 4X*  7HC0SS  =  *E15.8) 

c  check  sign  of  interproduct 

IF  (PR04 (J) )  442446 

4  LL  =  1 

14  J  =  J  +  1 

A( J)  =  A ( J— 1 )*AAA 
DO  Bis  24KCOUN 

5  U ( 1 4  I )  =  U( 2  4  I )  +  A ( J ) *S ( I ) 

PE TURN 
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6  I F (LL  • FO •  -1}  GO  TO  9 

C  CUBIC  INTFPPOL  AT'7  “OR  ALr’HA 

7  2  =  3*0* (COSS( U-l  )-COSS(U)  )/(A(  J)-A(U--1  )  >+c ‘•■vVi  (u) 

1+PR04J J-l ) 

W  =  SORT  {  (  J  )  04  ( J~ 1  3  5 

A  A  =  A  ( J  >  -  (  PROA  f  J  )  7  )  v-  (  A  (  J  )  - ' ■>  J-  !  5  '  /  <  '*  4  { 5  - 

1PR04{Uj1>  +  2.*X> 

24  DO  8  I  =  2 « JCCOUM 

8  U{  1  »  I  )  =  U( 2  «  I  3  +  AA  *F(  !  5 
LBO  =  -1 

NNZ»  =  20 
WRITE  (6.807)  A A 

802  FORMAT  <  2X «  LHAA  -  .  F 1  5  •  “  «  //  > 

‘ RETURN 

C  IF  FIRST  GUFSS  ON  At.*5 HA  GIVFS  POSITIVE 
C  G*S«  THEN  RFDUCF  ALPHA  GUFFS 

9  A  (  J)  =  #]  *A(J) 

DO  11  I  =  ?«KC OUN 

11  U(1  «  I  )  =  U< 2  .  I  )  + 

RETURN 

END 

S1F3FTC  SUFR 

SUBROUTINE  0FR1  (Hd) 

C0MV.0N/GC5CD/X(5*2001  3  »XX(3«">  «R''L  «F3  .DH<“>  , 
lG(2.?or>l  )  *  S  C  200 1  )  «U(2«2P0i  )  .Ft  7001  )  .  xr.Ouv  ♦  RRAV , 
2C0ST . N.N2P « L.RO •  M •  P » 0 1  . 07  i n't  .  o*'  ..v?.n 

DXtl « I )  r  XX (2.1) 

DX(2«n  =  ( 1  .+P)*0RAV*U<  :  »'•')  -  03*XX(1.1)  - 
1G4*XX(2.I)  -  08* ( r ( M ) —XX (3.1)1 
DX (3.1)  =  XX (4  « I ) 

DX ( 4 ,  J)  a  -R*0R6V*U ( 1  ♦ V }  +  06*  XX {1.15  +  C7*XX (2.1) 
1+05*{F<M)-XX(3. I > )  “  GRAV 
DX(5,I)  s"  (UM.lt)  -  Q1*XX<1.1)  -  C.?*XX  (2.15  -  1*0  ) 
1**2  +  P*U.(  1  «M)**2 
RETUPN 
END  . 

SIBFTC  SUBC 

SUBROUT  I NE  DFR2  ( H .  1) 

COMMON/GOCD/X ( 5 » 200 1 ) t  XX ( 5 « 5 ) .DX ( 5  «  3  3  »  DH ( 5 ) . 

1 G ( 2 . 200 1 ) .S(2001 ) ,U(2»?001 3 ."(2001 3 . <COU\ . G? AV . 
2C0ST « NN2P ♦ LBO  *M«  P . 0 1 ♦ 02  «  03  «  04 .OB  « 06 . 07 . o 
C  COMPUTE  GP An I ENT 

G(l.M)  =  (1 •+R)*GRAV*XX(2« I  3-  R*SPA V*XX ( 4 . I )  + 

1  2#*(U(1.M)  -  01*XC1»M)  -  02*X(?«M)  -  1.0) 

2+  2  •  -*P*U  (1  •  M ) 

DX(l.l)  =  03*XX(2«I>  -  OS*XX(4«I)  + 

1  2.*01*(U(  1  *M)  -  01*X(1,V.)  -  02*X(2«V)  -  -  1  •  C  3 
DX (2.1)  =.  -XX(l.I)  +  Q4*XX (2.1)  -  07* XX (4,1)  + 

1  2«*Q2*(U<  1  «M)  -  0 1 *X  ( 1  . M )  -  02*X(7«")  -  1**0) 

DX (3.1)  =  05* ( XX (4.1)  -  XX(2»i)) 

DX ( 4  « I )  =  —XX (3.  I ) 
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RETURN¬ 

ED 

£!oftc  ru^r 

SU3R0UT  I  \"~  Ri  JVG=-  |lJ*NilL) 

,CG!'-‘Vi1v/GrC.r)/Y<S«?Ori  )  «XX(?«r;  J  ♦DX(‘5«C>)  iDH(5)  » 
!G<5>»2001  )  •  S(?P»'l  1  .UC  5  «KCOUV,GP/-V* 

2C CRT « N:\-7P  « lp o •  m « °  •  r :  •  or-  •  or. « o* » of  » os « 07 « o 

C  4TH  ORPFR  RUNGE~<UTTA  !  N'TF  GOAT  I  CN 
N3  =  N  +  1 

NA  =  N  +  .7 

DH(N2)  =  H/2.0 

DH(N  +  2)=.  PH(N2> 

DPIN'M  =  H 
DO  20!  J  =  N7«  \'4 
DO  1 0!  !<  =  1  .  !L 

]0l  XX  (f»J)  =XX(  K  1  NS.)  +D-;  UVOXIKi  J— 1  > 

IF  <K>  .?«?*! 

1  CALL  PF-O 1  ( 11  ♦  J ) 

GO  TO  70 1 

2  CALL 

201  COvt  !vijF 

DO  AM  X  =  1  «  IL 

401  xx  { x  *n?  >  =xx  '  •* « n  >  -  t  t  < « v >  +? • -rx  c  x  «  n?.  >'+? .«r>x(K« n+p  ) 

1  -hDX  ( < « \'4 )  )/6>0 
IF  (H)  4,4,3 

3  V  =  M  -  1 

CALL  0*^1  <H«N?) 

GO  TO  F 

A  V.  -  M  -  J 

CALL  PFP2(K,\'?!. 

5  CON'T  1  \‘U~ 

RETURN 

END 

SIRLDR  FILES  IP-  DEC  69 

SFJLE  FILES  -UNIT07-«PP1 i READY «  OUTPUT, BCD* 3LK= 14 

SFDICT  FILES 

*4  P7G02-  *(-PPPPPP 

STEXT  FILES 

*4  =  4H*  I 

SCDICT  FILES 

*5  *M  3  1-*>P  3  XD- 

FDKEND  files 

T-EOF 
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